We  observe  that  solutions  of  a  large  class  of  highly  oscillatory  second  order  linear 
ordinary  differential  equations  can  be  approximated  using  nonoscillatory  phase 
functions.  In  addition,  we  describe  numerical  experiments  which  illustrate  com¬ 
putational  implications  of  this  fact.  For  example,  many  special  functions  of  interest 
—  such  as  the  Bessel  functions  Jv  and  Yv  —  can  be  evaluated  accurately  using  a 
number  of  operations  which  is  0(1)  in  the  order  v.  The  present  paper  is  devoted  to 
the  development  of  an  analytical  apparatus.  Numerical  aspects  of  this  work  will 
be  reported  at  a  later  date. 
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1.  Introduction 


Given  a  differential  equation 

y"(t)  +  A 2q(t)y(t)  =  0  for  all  0  ^  t  ^  1,  (1) 

where  A  is  a  real  number  and  q  :  [0, 1]  — *  M  is  smooth  and  strictly  positive,  a  sufficiently 
smooth  a  :  [0,1]  — >  M  is  a  phase  function  for  (1)  if  the  pair  of  functions  u,  v  defined  by  the 
formulas 


and 


u(t )  = 

cos  (a(f)) 

(2) 

1  a'«r 

v(t)  = 

sin(o:(t)) 

(3) 

|a'(i)|1/2 

form  a  basis  in  the  space  of  solutions  of  (1).  Phase  functions  have  been  extensively  studied: 
they  were  first  introduced  in  [9],  play  a  key  role  in  the  theory  of  global  transformations  of 
ordinary  differential  equations  [3, 10],  and  are  an  important  element  in  the  theory  of  special 
functions  [16,  6, 11, 1]. 


Despite  this  long  history,  an  important  property  of  phase  functions  appears  to  have  been 
overlooked.  Specifically,  that  when  the  function  q  is  nonoscillatory,  solutions  of  the  equation 
(1)  can  be  accurately  represented  using  a  nonoscillatory  phase  function. 


This  is  somewhat  surprising  since  a  is  a  phase  function  for  (1)  if  and  only  if  it  satisfies  the 
third  order  nonlinear  ordinary  differential  equation 


(«'W)2 


Wit) 


1  a'" {t)  3  fa"{t)  V 

2  a'{t)  +  4  ycPtf)  / 


for  all  0  <  t  <  1. 


(4) 


The  equation  (4)  was  introduced  in  [9],  and  and  we  will  refer  to  it  as  Rummer's  equation. 
The  form  of  (4)  and  the  appearance  of  A  in  it  suggests  that  its  solutions  will  be  oscillatory  — 
and  most  of  them  are.  However,  Bessel's  equation 

y{t)  =  0  for  all  0  <  t  <  oo  (5) 


y"{t)+  l- 


t2 


furnishes  a  nontrivial  example  of  an  equation  which  admits  a  nonoscillatory  phase  function 
regardless  of  the  value  of  A.  If  we  define  u,  v  by  the  formulas 


and 


(6) 

(7) 


where  J\  and  Y\  denote  the  Bessel  functions  of  the  first  and  second  kinds  of  order  A,  and  let 
a  be  defined  by  the  relations  (2),(3),  then 


a'(t)  = 


(8) 


nt  Jx(t)  +  Yx2(t)' 

It  can  be  easily  verified  that  (8)  is  a  strictly  negative,  monotonically  decreasing  function  of  t 
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on  (0,  oo).  That  this  nonoscillatory  phase  function  for  Bessel's  equation  exists  is  the  basis  of 
several  methods  for  the  evaluation  of  Bessel  functions  of  large  orders  and  for  the  computation 
of  their  zeros  [6,  8, 15]. 


The  general  situation  is  not  quite  so  favorable:  there  need  not  exist  a  nonoscillatory  function  a 
such  that  (2)  and  (3)  are  exact  solutions  of  (1).  However,  assuming  that  q  is  nonoscillatory  and 
A  is  sufficiently  large,  there  exists  a  nonoscillatory  function  a  such  that  (2),  (3)  approximate 
solutions  of  (1)  with  spectral  accuracy  (i.e.,  the  approximation  errors  decay  exponentially 
with  A). 


To  see  that  this  claim  is  plausible,  we  apply  Newton's  method  for  the  solution  of  nonlinear 
equations  to  Rummer's  equation  (4).  In  doing  so,  it  will  be  convenient  to  move  the  setting  of 
our  analysis  from  the  interval  [0, 1]  to  the  real  line  so  that  we  can  use  the  Fourier  transform 
to  quantity  the  notion  of  "nonoscillatory"  Suppose  that  the  extension  of  q  to  the  real  line  is 
smooth  and  strictly  positive,  and  that  the  Fourier  transform  of  log(g)  is  smooth  and  decays 
rapidly  Letting 

(a'(t))2  =  A2exp(r(f))  (9) 

in  (4)  yields  the  logarithm  form  of  Rummer's  equation: 

r"(t)  —  ^  (r'(t))2  +  4A2  (exp(r(f))  —  q{t ))  =  0  for  all  t  e  M.  (10) 

We  use  {rn}  to  denote  the  sequence  of  Newton  iterates  for  the  equation  (10)  obtained  from 
the  initial  guess 

r0(t)  =  log  (q(t)).  (11) 


The  function  r0  corresponds  to  the  first  order  WRB  approximations  for  (1).  That  is  to  say  that 
if  we  insert  the  associated  phase  function 

f  t  / 1  \  rt 

exp 


’t  / 1  \  rt 

ao(t )  =  A  exp  (  -r0(w)  )  du  =  A  \J q{u)du 

Jo  \2  /  Jo 


into  (2),(3),  then 


and 


u(t)  -  cos  ^A  J  \J q{u )  dv^j 


v(t )  =  q  1/4  (f)  sin  A 


a/^M 

Jo 


u )  du 


For  each  n  ^  0,  rn+ 1  is  obtained  from  rn  by  solving  the  linearized  equation 


where 


and  letting 


h"(t )  -  ~r'n(t)h'(t)  +  4A2  exp  (rn(t))  h(t)  =  fn(t )  for  all  t  e 


fn(t)  =  -r'n(t)  +  |  (r'Jt)f  -  4A2  (exp  (rn(t))  -  q(t )) , 


rn+i(t)  =  rn(t )  +  h(t). 


(12) 

(13) 

(14) 

(15) 

(16) 
(17) 
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By  introducing  the  change  of  variables 


"•t 

x(t)  = 

Jo 


exp 


rn(u) 


du 


(18) 

(19) 


into  (15),  we  transform  it  into  the  inhomogeneous  Helmholtz  equation 

h"(x)  +  4A 2h(x)  =  gn(x)  for  all 

where 

gn(x)  =  exp  (~rn(x))  fn(x).  (20) 

Suppose  that  gn  decays  rapidly  (when  n  =  0,  this  is  a  consequence  of  our  assumption  that 
log(g)  has  a  rapidly  decaying  Fourier  transform)  and  let  h*  be  the  solution  of  (19)  whose 
Fourier  transform  is 


MO  = 


MO 


(21) 


4A2  -  £2 

Since  h*(£)  is  singular  when  £  =  +2A,  h*  will  necessarily  have  a  component  which  oscillates 
at  frequency  2A.  However,  according  to  (21),  the  Lx  (M)  norm  of  that  component  is 

<7n(2A) 


In  fact,  by  rearranging  (21)  as 


MO  =  7T 


4A 


i  /MO  ,  MO 


+ 


4A  \2A  —  £  2A  +  £ 

and  decomposing  each  of  the  terms  on  the  right-hand  side  of  (23)  as 

MO  1 


2A  +  £  4A 

we  obtain 


MO  -^.(+2A)exp(-(2A  +  02)  ,  -  ^oUexp  (-(2A  ±  00 


2A  +  £ 

h*(x)  =  h0(x )  +  hi(x), 


2A  +  £ 


where  h0  is  defined  by  the  formula 


1  ( 9n(0  -  9n(-2X)  ex p(-(2A  +  02)  ,  MO  “  M 2A)exp(-(2A-£)2) 


h°^  4A  V  2 A  +  £ 

and  hi  is  defined  by  the  formula 


+ 


MO  —  (  9n{  2A) 


exp(-(2A  +  02) 


+ .(/» (2  A) 


2A^£ 

exp(-(2A^02) 


(22) 

(23) 

(24) 

(25) 

(26) 

(27) 


2A  +  £  v  '  2A  —  £ 

Since  the  factor  in  the  denominator  in  (26)  has  been  canceled  and  both  gn  and  the  Gaussian 
function  are  smooth  and  rapidly  decaying,  h0  is  also  smooth  and  rapidly  decaying.  Mean¬ 
while,  a  straightforward  calculation  shows  that  the  Fourier  transform  of 

4  6lf  (f  )  exP(2AM  (28) 


is 


exp(-(2A-Q2) 

2A-£ 


(29) 
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(30) 


so  that  (27)  implies  that 

hi(x)  =  ^(-2A)^r  erf  exp(2A ix)  -  .gn(2A)^erf'  exp(-2Aix) 

Since  r/n  is  real-valued,  gn{ 2A)  =  ( — 2 A ) .  Inserting  this  into  (30)  yields 

erf  sin  (2Ax) ,  (31) 

which  makes  it  clear  that  the  L00  (M)  norm  of  hi  is  (4A)_1g),,(2A). 

In  (25),  the  solution  of  (19)  is  decomposed  as  the  sum  of  a  nonoscillatory  function  h0  and  a 
highly  oscillatory  function  h  \  of  small  magnitude.  However,  the  solution  of  (15)  is  actually 
given  by  the  function 

h*(x(t ))  =  ho(x(t))  +  hi(x(t))  (32) 

obtained  by  reversing  the  change  of  variables  (18).  But  since  x(t)  is  nonoscillatory  and  the 
Fourier  transform  of  h0(x)  decays  rapidly,  we  expect  that  the  composition  h0(x(t ))  will  also 
have  a  rapidly  decaying  Fourier  transform.  The  L00  (M)  norm  of  hi(x(t))  is,  of  course,  the 
same  as  that  of  h\  (. x ).  So  the  solution  of  the  linearized  equation  (15)  can  be  written  as  the  sum 
of  a  nonoscillatory  function  h0(x(t ))  and  a  highly  oscillatory  function  hi  (x(t))  of  negligible 
magnitude. 

If,  in  each  iteration  of  the  Newton  procedure,  we  approximate  the  solution  of  (15)  by  con¬ 
structing  h*(x(t ))  and  discarding  the  oscillatory  term  hi(t(x))  of  small  magnitude,  then  it  is 
plausible  that  we  will  arrive  at  an  approximate  solution  r(t)  of  the  logarithm  form  of  Rum¬ 
mer's  equation  which  is  nonoscillatory,  assuming  the  Fourier  transform  of  r0(f)  =  log (q(t)) 
decays  rapidly  enough  and  A  is  sufficiently  large. 

Most  of  the  remainder  of  this  paper  is  devoted  to  developing  a  rigorous  argument  to  replace 
the  preceding  heuristic  discussion.  In  Section  2,  we  summarize  a  number  of  well-known 
mathematical  facts  to  be  used  throughout  this  article.  In  Section  3,  we  reformulate  Rummer's 
equation  as  a  nonlinear  integral  equation.  Once  that  is  accomplished,  we  are  in  a  position  to 
state  the  principal  result  of  the  paper  and  discuss  its  implications;  this  is  done  in  Section  4. 
The  proof  of  this  principal  result  is  contained  in  Sections  5,  6,  7  and  8.  Section  9  contains  an 
elementary  proof  of  a  relevant  error  estimate  for  second  order  differential  equations  of  the 
form  (1). 

In  Section  10,  we  present  the  results  of  numerical  experiments  concerning  the  evaluation  of 
special  functions.  The  details  of  our  numerical  algorithm  will  be  reported  at  a  later  date. 

We  conclude  with  a  few  brief  remarks  in  Section  11. 


2.  Preliminaries 

2.1.  Modified  Bessel  functions 

The  modified  Bessel  function  Kv{t)  of  the  first  kind  of  order  v  is  defined  for  lei  and  v  e  C 
by  the  formula 

r*CO 

Kv(t)  =  exp  (— f  cosh  (f))  cosh(z/f)  dt.  (33) 

Jo 


5 


The  following  bound  on  the  ratio  of  Ku+L  to  K„  can  be  found  in  [14]. 


Theorem  1.  Suppose  that  t  >  0  and  is  >  0  are  real  numbers.  Then 

Ku+i(t)  is  +  V  is2  +  f 2  2  is 

Kv{t)  <  t  t  + 


(34) 


2.2.  T/?e  binomial  theorem 

A  proof  of  the  following  can  be  found  in  [13],  as  well  as  many  other  sources. 

Theorem  2.  Suppose  that  r  is  a  real  number,  and  that  y  is  a  real  number  such  that  \y\  <  1.  Then 


(1+9)'=  £ 


T(r  +  1) 


-f  T(k  +  1  )T(r  -  k  +  1) 


(35) 


k= 0 


2.3.  The  Lambert  W  function 

The  Lambert  W  function  or  product  logarithm  is  the  multiple-valued  inverse  of  the  function 

f(z)  =  zexp(z).  (36) 

We  follow  [5]  in  using  Wo  to  denote  the  branch  of  W  which  is  real-valued  and  greater  than 
or  equal  to  —1  on  the  interval  [—1/e,  —  oo)  and  W-±  to  refer  to  the  branch  which  is  real-valued 
and  less  than  or  equal  to  —1  on  [—1/e,  0). 

We  will  make  use  of  the  following  elementary  facts  concerning  W0  and  W_i  (all  of  which  can 
be  found  in  [5]  and  its  references). 

Theorem  3.  Suppose  that  y  ^  —1/e  is  a  real  number.  Then 

rrexp  (x)  ^  y  if  and  only  if  x^Wo(y).  (37) 

Theorem  4.  Suppose  that  0  <  y  ^  1/e  is  a  real  number.  Then 

xexp(— x)  ^  y  if  and  only  if  x  ^  —W-i(-y).  (38) 


Theorem  5.  Suppose  that  0  ^  x  ^  1  is  a  real  number.  Then 

|  <  Wq(x)  ^  X. 


Theorem  6.  Suppose  that  x  >  exp(l)  is  a  real  number.  Then 


log(a;)  ^  —  W-i 


^  2  log(x). 


(39) 


(40) 
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2.4.  Frechet  derivatives  and  the  contraction  mapping  principle 


Given  Banach  spaces  X,  Y  and  a  mapping  /  :  X  — »  Y  between  them,  we  say  that  /  is  Frechet 
differentiable  at  x  e  A"  if  there  exists  a  bounded  linear  operator  X  — >  V’,  denoted  by  /',  such 
that 


\\f(x  +  h)-f(x)-fx[h]\\ 
h->  o  ||  fi|| 


0. 


(41) 


Theorem  7.  Suppose  that  X  and  Y  are  a  Banach  spaces  and  that  f  :  X  —*  Y  is  Frechet  differentiable 
at  every  point  of  X.  Suppose  also  that  D  is  a  convex  subset  of  X,  and  that  there  exists  a  real  number 
M  >  0  such  that 


for  all  x  e  D.  Then 


for  all  x  and  y  in  D. 


\\a  <  m 

f(x)  —  f(y)\\  <  Af ||a:  —  j/|| 


(42) 

(43) 


Suppose  that  f  :  X  — *  X  is  a  mapping  of  the  Banach  space  X  into  itself.  We  say  that  /  is 
contractive  on  a  subset  D  of  X  if  there  exists  a  real  number  0  <  a  <  1  such  that 


\\f(x)  ~  f(y)\\  <  a\\x-y\\  (44) 

for  all  x,y  e  D.  Moreover,  we  say  that  {xn}f=0  is  a  sequence  of  fixed  point  iterates  for  /  if 

Xn+1  =  f(xn)  for  all  n  ^  0. 


Theorem  7  is  often  used  to  show  that  a  mapping  is  contractive  so  that  the  following  result 
can  be  applied. 


Theorem  8.  (The  Contraction  Mapping  Principle)  Suppose  that  D  is  a  closed  subset  of  a  Banach 
space  X.  Suppose  also  that  f  :  X  — >  X  is  contractive  on  D  and  f(D)  c  D.  Then  the  equation 

x  =  f(x)  (45) 

has  a  unique  solution  a*  e  D.  Moreover,  any  sequence  of  fixed  point  iterates  for  the  function  f  which 
contains  an  element  in  D  converges  to  a*. 

A  discussion  of  Frechet  derivatives  and  proofs  of  Theorems  7  and  8  can  be  found,  for  instance, 
in  [18]. 


2.5.  Gronwall's  inequality 

The  following  well-known  inequality  can  be  found  in,  for  instance,  [2]. 

Theorem  9.  Suppose  that  f  and  g  are  continuous  functions  on  the  interval  [a,  b]  such  that 

f(t )  0  and  g(t)  ^  0  for  all  a  <  £  <6.  (46) 

Suppose  further  that  there  exists  a  real  number  C  >  0  such  that 

f(t )  <  C  +  f  f(s)g(s)  ds  for  all  a  <  £  <6.  (47) 
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Then 


fit)  Cexp 


g{s)  ds 


for  all  a  ^  t  ^  b. 


(48) 


2.6.  Schwarzian  derivatives 


The  Schwarzian  derivative  of  a  smooth  function  /  :  M  — *  R.  is 

{/ 1}  _  n*>  _  3  rm \2  (49) 

u' !  /'to  2  v  m)  ' 

If  the  function  x(t)  is  a  diffeomorphism  of  the  real  line  (that  is,  a  smooth,  invertible  mapping 
M  — >  M),  then  the  Schwarzian  derivative  of  x(t)  can  be  related  to  the  Schwarzian  derivative  of 
its  inverse  t(x);  in  particular. 


The  identity  (50)  can  be  found,  for  instance,  in  Section  1.13  of  [11]. 


(50) 


3.  Integral  equation  formulation 

In  this  section,  we  reformulate  Rummer's  equation 

,  ..  ..2  o  ,  \  1  a'"(t)  3  / a" (t)\2 

(a(t»  2zw+ *  m) 

as  a  nonlinear  integral  equation.  As  in  the  introduction,  we  assume  that  the  function  q  has 
been  extended  to  the  real  line  and  we  seek  a  function  a  which  satisfies  (51)  on  the  real  line. 


(51) 


By  letting 


in  (51),  we  obtain  the  equation 


(n'(f))2  =  A2exp(r(t)) 


r"(t)  —  -  ( r'(t ))  +  4 A2  (exp(r(i))  —  q(t))  =  0  for  all  t  e 

We  next  take  r  to  be  of  the  form 

r(t)  =  log  (q(t))  +  S(t), 

which  results  in 

d"(t)  -  -  \  id'it))2  +  4A 2q(t)  (exp (6(f))  -  1)  =  q(t)p(t),  fo 

where  p  is  defined  by  the  formula 

1  /  5  /  q'(t)\2  q"(t ) 


p{t)  = 


q(t)  V 4  V  q(t) 


q(t) 


(52) 

(53) 

(54) 

(55) 

(56) 


Note  that  the  function  p  appears  in  the  standard  error  analysis  of  WKB  approximations  (see, 
for  instance,  [12]).  Expanding  the  exponential  in  a  power  series  and  rearranging  terms  yields 
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the  equation 

s'\t) - + 4x^mt)  - 1  mf + +  •••')  =  ««?(*)•  (57) 


Applying  the  change  of  variables 

x(t)  =  f  v 

Jo 

/q(u)  du 

(58) 

transforms  (57)  into 

5'\x)  +  4A 2d(x)  -  \  L(M  +  4A2 

A(A  + 

{KM 

3! 

'j  =  p{x)  for  all  x  eR.  (59) 

At  first  glance,  the  relationship  between  the  function  p(x)  appearing  in  (59)  and  the  coefficient 
q(t)  in  the  ordinary  differential  equation  (1)  is  complex.  However,  the  function  pit)  defined 
via  (56)  is  related  to  the  Schwarzian  derivative  (see  Section  2.6)  of  the  function  x(t)  defined  in 
(58)  via  the  formula 

It  follows  from  (60)  and  Formula  (50)  in  Section  2.6  that 

p(x)  =  2{t,x}.  (61) 

That  is  to  say:  p,  when  viewed  as  a  function  of  x,  is  simply  twice  the  Schwarzian  derivative 
of  t  with  respect  to  x. 

It  is  also  notable  that  the  part  of  (59)  which  is  linear  in  S  is  a  constant  coefficient  Helmholtz 
equation.  This  suggests  that  we  form  an  integral  equation  for  (59)  using  a  Green's  function 
for  the  Helmholtz  equation.  To  that  end,  we  define  the  integral  operator  T  for  functions 
/  e  L1  (M)  via  the  formula 

r[/](x)  =  ^J  sin  (2  \  \x -y\)  f(y)  dy  (62) 

The  following  theorem  summarizes  the  relevant  properties  of  the  operator  T. 

Theorem  10.  Suppose  that  A  >  0  is  a  real  number,  and  that  the  operator  T  is  defined  as  in  (62). 
Suppose  also  that  f  e  L1  (M)  n  C  (M).  Then: 

1.  T  [/]  (x)  is  an  element  of  C2  (M); 

2.  T  [f]  (x)  is  a  solution  of  the  ordinary  differential  equation 

y"(x)  +  4A  2y(x)  =  f(x)  for  all  ieR;  and 

3.  the  Fourier  transform  ofT[f]  ( x )  is  the  principal  value  of 

m  =±( ML  .  ML) 

4A2-e  4A  \^2A  —  £  2A  +  f)  ‘ 
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Proof.  We  observe  that 

rmw-i 


sin  (2A  (x  —  y))  f{y)  dy  +  — 


=  —  sin  (2  Ax 
4A  v 


V-  - 


—  00 
00 


cos(2A y)f(y)  dy  -  —  cos(2Ax 


sin  (2X(y  -  x))  f(y)  dy 

sin(2A y)f(y)  dy  (63) 


—00 

'00 


+  —  cos 
4A 


•»oo  -|  ro o 

(2Ax)  sin(2A y)f(y)  dy  -  —  sin(2Ax)  cos (2 \y)  f(y)  dy 

Jx  Jx 


for  all  x  e  M.  We  differentiate  both  sides  of  (63)  with  respect  to  x,  apply  the  Lebesgue  dom¬ 
inated  convergence  theorem  to  each  integral  (this  is  permissible  since  the  sine  and  cosine 
functions  are  bounded  and  /  e  L1  (M))  and  combine  terms  in  order  to  conclude  that  T  [/]  is 
differentiable  everywhere  and 


d  1  r°° 

— 2 


cos  (2A  \x  -  y\)  sign(x  -  y)f(y)  dy 


for  all  x  s  M.  In  the  same  fashion,  we  conclude  that 

d 


-)  T  [f  ]  (x)  =  f(x)  —  X 


sin  (2A  \x  -  y\)  f(y )  dy 


(64) 


(65) 


for  all  x  e  M.  Since  /  is  continuous  by  assumption  and  the  second  term  appearing  on  the 
right-hand  side  in  (65)  is  a  continuous  function  of  x  by  the  Lebesgue  dominated  convergence 
theorem,  we  see  from  (65)  that  T  [  f  ]  is  twice  continuously  differentiable.  By  combining  (65) 
and  (62),  we  conclude  that  T  [f]  is  a  solution  of  the  ordinary  differential  equation 

y"(x)  +  4A 2y(x)  =  f{x)  for  all  x  e  M.  (66) 


We  now  define  the  function  g  through  the  formula 


9(0  =  77 


+ 


(67) 


4A  \2A  —  £  2A  +  f, 

It  is  well  known  that  the  Fourier  transform  of  the  principal  value  of  1/x;  is  the  function 

— i7rsign(x);  (68) 

see,  for  instance,  [17]  or  [7].  It  follows  readily  that  the  inverse  Fourier  transform  of  the  prin¬ 
cipal  value  of 

1  (69) 

(70) 


2A  +  £ 


is 


+  —  exp  (+2A ix)  sign(x). 

j£j  i 


From  this  and  (67),  we  conclude  that 


1/1  1 

g(x)  =  —  I  —  exp  (— 2A ix)  sign(x)  —  —  exp  (2A ix)  sign(x) 

Tt/\  \  £  i  £  i 

=  ^sin(2A|x|). 


(71) 
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In  particular,  T  [/]  is  the  convolution  of  /  with  g.  As  a  consequence, 

1 


T[fm  =  9(0f(0  = 


m  m 


4A  y2A  —  £  2A  +  £ 

which  is  the  third  and  final  conclusion  of  the  theorem. 

In  light  of  Theorem  10,  it  is  clear  that  introducing  the  representation 

S(x)  =  T[a]  (x) 

into  (59)  yields  the  nonlinear  integral  equation 

a(x)  =  S  [T  [cr]  ]  (x)  +  p(x)  for  all  x  e  M, 
where  S  is  the  operator  defined  for  functions  /  e  C1  (M)  by  the  formula 


4A 


2 ,  (/Mr  +  Ml  +  [Ml  + 


2! 


3! 


4! 


(72) 

□ 

(73) 

(74) 

(75) 


The  following  theorem  is  immediately  apparent  from  the  procedure  used  to  transform  Rum¬ 
mer's  equation  (51)  into  the  nonlinear  integral  equation  (74). 

Theorem  11.  Suppose  that  A  >  0  is  a  real  number,  that  q  :  M  — »  M  is  an  infinitely  differentiable, 
strictly  positive  function,  that  x(t)  is  defined  by  (58),  and  that  p(x)  is  defined  via  (61).  Suppose  also 
that  a  e  L1  (M)  n  C  (M)  is  a  solution  of  the  integral  equation  (74),  that  5  is  defined  via  the  formula 

1  f00 

S(x)  =  T[a]  (x)  =  —  J  sin(2A|x  -y\)a(y)  dy,  (76) 

and  that  the  function  a  is  defined  by  the  formula 

aft)  =  A  J  ^/q(u)exp  du.  (77) 


Then: 

1.  5(x )  is  a  twice  continuously  differentiable  solution  of  (59); 

2.  5(x(t))  is  a  twice  continuously  differentiable  solution  of  of  (57); 

3.  a(t)  is  three  times  continuously  differentiable  solution  of  (51);  and 

4.  a(t)  is  a  phase  function  for  the  ordinary  differential  equation 

y"(t)  +  A 2q(t)y(t)  =  0  for  all  0  <  f  <  1.  (78) 


4.  Overview  and  statement  of  the  principal  result 

The  composition  operator  S  o  T  appearing  in  (74)  does  not  map  any  Lebesgue  space  L"  (M)  or 
Holder  space  Cfc,Q(M)  to  itself,  which  complicates  the  analysis  of  (74).  Moreover,  the  integral 
defining  T  [u]  only  exists  if  either  a  e  L1  (M)  or  u(+2A)  =  0.  Even  if  both  of  these  conditions 
are  satisfied,  it  is  not  necessarily  the  case  that  S  [T  [cr]  ]  +  p  will  satisfy  either  condition.  This 
casts  doubts  on  whether  (74)  is  solvable  for  arbitrary  p. 
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We  avoid  a  detailed  discussion  of  which  spaces  and  in  what  sense  (74)  admits  solutions  and 
instead  show  that,  under  mild  conditions  on  p  and  A,  there  exists  a  function  pb  which  approx¬ 
imates  p  and  such  that  the  equation 

a(x)  =  S[T[a]](x) +pb(x)  for  all  x  e  M  (79) 

admits  a  solution.  Moreover,  we  prove  that  if  p  is  nonoscillatory  then  the  solution  of  (79)  is 
also  nonoscillatory  and  \\p  —  Pb\\ao  decays  exponentially  in  A.  The  next  theorem,  which  is  the 
principal  result  of  this  article,  makes  these  statements  precise.  Its  proof  is  given  in  Sections  5, 
6,  7  and  8. 


Theorem  12.  Suppose  that  q  e  C 00  (M)  is  a  strictly  positive,  and  that  x(t)  is  defined  by  the  formula 


-*£ 

x(t)  =  Vol 

Jo 


u)  du. 


(80) 


Suppose  also  that  p(x)  is  defined  via  the  formula 

p(x)  =  2{t,x}]  (81) 

that  is,  p(x)  is  twice  the  Schwarzian  derivative  of  the  variable  t  with  respect  to  the  variable  x  defined 
via  (80).  Suppose  furthermore  that  there  exist  real  numbers  A  >  0,  T  >  0  and  a  >  0  such  that 


A  ^  6  max  <  — ,  T 
a 


(82) 


and 


^  Texp  (^a  |£|)  for  all  (el.  (83) 

Then  there  exist  functions  pb  e  C00  (M)  and  ab  e  L'2  (M)  n  C*00  (M)  such  that  ab  is  a  solution  of  (79), 

|ab(0|  <  2T  exp  for  all  £  e  R,  (84) 


and 


24r 


P  —  P&||ao  <  —  exP  -paX  . 


5  a 


5 


6 


(85) 


According  to  Theorem  11,  if  a  is  a  solution  of  the  integral  equation  (74),  then  the  function  a 
defined  by  (77)  is  a  phase  function  for  the  differential  equation  (1).  We  define  ab  in  analogy 
with  (77)  using  the  solution  ab  of  the  modified  equation  (79).  That  is,  we  let  5b  be  the  function 
defined  by  the  formula 

1  f°° 

5b(x)  =  T  [cr6]  (x)  =  —  J  sin(2A|x  -  y\)ab(y)  dy ,  (86) 


and  then  define  ab  via 

(X-bif)  =  A  J  VVMexP  du-  (87) 

Obviously,  ab  is  not  a  phase  function  for  the  equation  (1).  However,  if  we  define  qb  :  R  —>  R 
by  the  formula 


J_(5  fq_lftf\2 

Qb(t)  \4:  \qb(t)  J 


Pb(t), 


(88) 
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then  Theorem  11  implies  that  ab  is  a  phase  function  for  the  differential  equation 

y"(t)  +  A 2qb(t)y(t)  =  0  for  all  0  <  £  <  1.  (89) 

Since  \\q  —  gb||oo  is  bounded  in  terms  of  || p  —  pb\\ Xl  and  the  solutions  of  (89)  closely  approximate 
those  of  (1)  when  ||g  —  gb||oo  is  small,  the  function  ab  can  be  used  to  approximation  solutions 
of  (1)  when  || p  —  pfe||oo  is  small. 

Theorem  13,  which  appears  below,  gives  a  relevant  error  estimate.  Given  e  >  0,  it  specifies  a 
bound  on  || p  —  pb\\x  which  ensures  that  solutions  of  (89)  approximate  those  of  (1)  to  relative 
precision  e.  Its  proof  appears  in  Section  9. 


Defintion  1.  We  say  that  a  is  an  e-approximate  phase  function  for  the  ordinary  differential  equation 
(1)  if  there  exists  a  basis  of  solutions  {u,  v}  of  (1)  such  that 

| u(t)  —  u(t) |  ^  e  sup  \u(t)\  for  all  0  <  f  ^  1  (90) 


and 


v(t) 


where  u,  v  are  defined  by 


and 


v(t) |  ^  e  sup  \v(t)\  for  all  0  ^  t  ^  1, 
o=st^i 


u{t) 


sin(o;(f)) 


v{t) 


cos(a(£)) 

l«'(t)|1/2‘ 


(91) 


(92) 


(93) 


Theorem  13.  Suppose  that  q  e  G03  (M)  is  strictly  positive,  that  the  function  p  is  defined  by  the 
formula  (56),  and  that  there  exist  real  numbers  0  <  ip  <  y2  such  that 

Tj\  ^  q(t)  ^  rj2  for  all  0  ^  t  ^  1, 

W(t)\ 


q(t) 


^  rl'2  for  all  0  <  t  <  1, 


and 


\p(t)\  ^  rj2  for  all  0  <  £  <  1. 
Suppose  also  that  A  >  0,  e  >  0  are  real  numbers  such  that 


0  <  e  <  A  exp 

and  that  k  is  the  real  number  defined  by 


3/4' 

vf 


V2 


Jh 


Suppose  furthermore  that  pb  : 


k  =  20  —  +8  rji  +  10—  +  1. 

\Vi )  Vi 

— >  M  is  an  infinitely  differentiable  function  such  that 


V  ~  PbU  <  ~  exp  (~k)  exp  (  - 


3/4' 

W 


(94) 

(95) 


(96) 


(97) 


(98) 


(99) 
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Then  there  exists  a  function  qb  :  [0, 1]  — >•  M  such  that  (88)  holds  and  any  phase  function  for  (89)  is  an 
e-approximate  phase  function  for  (1). 

By  combining  Theorems  11, 12  and  13  we  obtain  the  following. 


Corollary  1.  Suppose  that  q,  x(t),  p,  ry,  q2,  A,  e,  k,  T,  a  and  q  are  as  in  the  hypotheses  of  Theorem  12 
and  13.  Suppose  that,  in  addition, 


A  >i-\k  + 
5  a 


3/4 

V2 


+  log 


/mM 

\?7i az  e 


Then  there  exists  a  function  ab  e  L 2  (M)  n  C00  (M)  such  that 

I  of.  (01  <  2Texp  for  all  £ 

and  the  function  ab(t )  defined  by 


o;6(£)  =  aJ  x/q(u)  exp  da, 

where  8b  is  the  function  defined  through  the  formula 

Mx)  =  ^  J  sin(2A|x  -  y\)ab{y)  dy , 

is  an  e-approximate  phase  function  for  the  ordinary  differential  equation  (1). 


(100) 


(101) 


(102) 


(103) 


Remark  1.  The  hypothesis  (100)  can  be  replaced  with  the  weaker  condition 


A  >  f  W-!  (  - 
5a 


25  qi a2 


exp  (—k)  exp  — 


(104) 


where  4b_i  denotes  the  branch  of  the  Lambert  W  function  which  is  real-valued  and  less  than  —1  on 
the  interval  [—1/e,  0)  (see  Section  2.3). 


The  proof  of  Theorem  12  is  divided  amongst  Sections  5,  6,  7  and  8.  The  principal  difficulty 
lies  in  constructing  a  function  pb  which  approximates  p  and  for  which  (79)  admits  a  solution. 
We  accomplish  this  by  introducing  a  modified  integral  equation 

a(x)  =  S[Tb[a]](x)  +  p(x),  (105) 

where  Tb  is  a  "band-limited"  version  of  T.  That  is,  Tb  [  f  ]  is  defined  via  the  formula 

TIT]©  =  nTmm,  w 

where  &(£)  is  a  C'f  (M)  bump  function.  This  modified  integral  equation  is  introduced  in  Sec¬ 
tion  5. 


In  Section  6,  we  show  that  under  mild  conditions  on  p  and  A,  (105)  admits  a  solution  a. 
The  argument  proceeds  by  applying  the  Fourier  transform  to  (105)  and  using  the  contraction 
mapping  principle  to  show  that  the  resulting  equation  admits  a  solution. 

In  Section  7,  we  estimate  the  Fourier  transform  of  the  solution  a  of  (105)  under  the  assumption 
that  p  is  exponentially  decaying.  We  show  that  a  is  also  exponentially  decaying,  albeit  at  a 
slightly  slower  rate. 
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In  Section  8,  we  define  a  band-limited  version  a b  of  the  solution  a  of  (105)  via  the  formula 

=  mm-  ao7) 

By  combining  the  observation  that  Tb  \  a  \  =  T  \  ab  |  with  (105),  we  obtain 

a(x)  =  S[T[ab]](x) +p(x).  (108) 

We  then  define  pb  by  the  formula 

pb{x)  =  p{x)  +  crb(x)  -  a(x)  (109) 


and  rearrange  (108)  as 

ab(x)  =  S[T[ab]]  (x)  +pb(x). 

The  decay  estimate  on  a  is  then  used  to  bound  \\p  —  pb\\oo  =  ||cr  —  crfe||oo  <  ||a  —  a), lx. 


(110) 


5.  Band-limited  integral  equation 

In  this  section,  we  introduce  a  "band-limited"  version  of  the  operator  T  and  use  it  to  form  an 
alternative  to  the  integral  equation  (74). 

Let  6(f)  be  any  infinitely  differentiable  function  such  that 

1.  6(f)  =  1  for  all  |f  |  s$  A, 

2.  0  ^  6(f)  ^  1  for  all  A  <  |f  |  <  v'2A,  and 

3.  6(f)  =  0  for  all  |f|  ^  V2X. 

We  define  Tb  [  f  ]  for  functions  /  e  L 1  (M)  via  the  formula 

Sm©  -  an) 

We  will  refer  to  Tb  as  the  band-limited  version  of  the  operator  T  and  and  we  call  the  nonlinear 
integral  equation 

a(x)  =  S[Tb[a]](x)  +  p(x)  for  all  x  e  M  (112) 

obtained  by  replacing  T  with  Tb  in  (74)  the  "band-limited"  version  of  (74). 

Since  Tb  is  a  Fourier  multiplier,  it  is  convenient  to  analyze  (112)  in  the  Fourier  domain  rather 
than  the  space  domain.  We  now  introduce  notation  which  will  allow  us  to  write  down  the 
equation  obtained  by  applying  the  Fourier  transform  to  both  sides  of  (112). 


We  let  Wb  and  Wb  be  the  linear  operators  defined  for  /  e  L 1  (M)  via  the  formulas 


m  [/]  (£)  -  rji 

(113) 

and 

(L  [/]«)  =  /«) 

(114) 

where  6(f)  is  the  function  used  to  define  the  operator  Tb. 
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For  functions  /  e  L1  (M),  it  is  standard  to  denote  the  Fourier  transform  of  the  function 

exp (f(x))  by  exp*  [/];  that  is, 

exp*  [/]  «)  =  ,5(0  +  m  +  hiM  +  l  *  R  /K)  +  ■  ■  ■  .  (115) 

In  (115),  <5  defers  to  the  delta  distribution  and  /*/*•••*/  denotes  repeated  convolution  of 
the  function  /  with  itself.  The  Fourier  transform  of  exp(/(x))  never  appears  in  this  paper; 
however,  we  will  encounter  the  Fourier  transforms  of  the  functions 

exp (f(x))  -  1  (116) 

and 

exP (f(x))  -  f(x)  -  1.  (117) 

So,  in  analogy  with  the  definition  (115),  we  define  exp/  [/]  for  /  el1  (M)  by  the  formula 

exp(  [/]  (?)  =  /( o  +  MO  + ;  —  +■■■•  (U8) 

and  we  define  exp/  [/]  for  /  e  L  l  (M)  via  the  formula 

exp;[/](0  =  ^  +  ^M  +  -'.  (119) 

That  is,  exp/  [/]  is  obtained  by  truncating  the  leading  term  of  exp*  [/]  and  exp/  [/]  is  obtained 
by  truncated  the  first  two  leading  terms  of  exp*  [/]. 

Finally,  we  define  functions  -0(£)  and  v(£)  using  the  formulas 

m  =  (120) 

and 

v(o=m-  (i2i) 

Applying  the  Fourier  transform  to  both  sides  of  (112)  results  in  the  nonlinear  equation 

m  =  R^  m  (122) 

where  R  [f]  is  defined  for  /  el1  (M)  by  the  formula 

R [/]«)  =  [/]  •  mu]  ©  - 4A2 exP; [wb [/] ]  (?)  + ««).  (123) 


6.  Existence  of  solutions  of  the  band-limited  equation. 

In  this  section,  we  give  conditions  under  which  the  sequence  {v'n}n=o  °f  fixed  point  iterates  for 
(122)  obtained  by  using  the  function  v  defined  by  (121)  as  an  initial  approximation  converges. 
More  explicitly,  T'0  is  defined  by  the  formula 

MO  =  v(0,  (124) 

and  for  each  integer  n  ^  0,  wn+  \  is  obtained  from  u>n  via 

Mi(0  =  R[M(0-  (125) 
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(126) 


Theorem  14.  Suppose  that  A  >  0  is  a  real  number,  and  that  v  e  L 1  (M)  such  that 

A2 

H  i  E  — . 
ll  in  lg 

Theen  the  sequence  {ipn)n=o  defined  by  (124)  and  (125)  converges  in  L 1  (M)  norm  to  a  function  ip 
which  satisfies  the  equation  (122)  for  almost  all  (el. 

Proof  We  observe  that  the  Frechet  derivative  (see  Section  2.4)  of  R  at  /  is  the  linear  operator 
R'f  :  L1  (E)  — >  L 1  (E)  given  by  the  formula 


R'f  [h]  (0  = 


Wb[f]*Wb[h]{(f) 


4A2  exp*  [Wb[f]\*Wb[h](0. 


From  formulas  (113)  and  (114)  and  the  definition  of  5(£)  we  see  that 


and 


immi^ 


Wb[f] 


\\fh 

2A2 

<  \\fh 

i  "  V2A 


(127) 

(128) 

(129) 


for  all  /  el1  (M).  From  (127),  (128)  and  (129)  we  conclude  that 

II . «'/  M II  i « j  | ' m  l f 1 1 1 1|  i w,,  [, ft  ]  ||  +  i. \2  n  ■ w„  [ f]\ || . ! exp  (j ||  ■ in [ . r ll  i : J )  ||  ■ m  i ft] | , , 


< 


< 


\\f\m 


4A2 


-  +  4A: 


:  l/|l  Ill'll 


l^+Mlexp 


2A2  2A2 

ll/lli 


exp 


I/ll 


4A2 


A2 


2A2 


2A2 

Hi 


(130) 


for  all  /  and  h  in  L 1  (M).  Similarly,  by  combining  (123),  (128)  and  (129)  we  conclude  that 

II R  [  /  ]  II  i «  \  |  w,  [  /  ]  [  +  W  ||  wb  [  /  ]  i?  exP  ( ||  m  [  /]  || , )  +  h  i 


g!ZE  +  M|expn/lli 


8A2  A2 


2A2 


(131) 


+  Mi 


whenever  /  el1 


Now  we  set  r  =  A2/9  and  let  B  denote  the  closed  ball  of  radius  r  centered  at  0  in  L 1 

l<wjl 


.  Since 
(132) 


where  Wo  denotes  the  branch  of  the  Lambert  W  function  which  is  real-valued  and  greater 
than  or  equal  to  —1  on  the  interval  [—1/e,  oo)  (see  Section  2.3),  we  conclude  that 

1 
4 

According  to  Theorem  5  of  Section  2.3  and  (132),  it  also  the  case  that 


r  /  r  \  1 

PeXPfej<- 


r  1 

A2  <  8" 


(133) 

(134) 
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^  1  —  +  —  exP  (  T7t))  Nil  <  Q2  +  i)  Nil  <  ^N 


(  r  r  (  r  y 

\w>  +  -’exp\w). 


We  combine  (133),  (134)  with  (130)  to  conclude  that 

KMII, 

for  all  h  e  L 1  (M)  and  /  e  B.  In  other  words,  the  L 1  (M)  — »  L 1 
operator  Rj  is  bounded  by  1/2  whenever  /  is  in  the  ball  B. 

Similarly,  we  insert  (133),  (134)  and  (126)  into  (131)  to  conclude  that 

2  2 

-  -I-  —  pyn  I  -  I  +  — 


(135) 

operator  norm  of  the  linear 


+  V2W 


<  r 


r  r  /  r  \ 

8P  +  A^eXP  M 


1 

+  2 


(136) 


111 

^  r  | - 1 - 1 — 

.64  4  2 

<  r 


for  all  /  e  B. 

Together  with  Theorem  7  of  Section  2.4,  formula  (135)  implies  that  the  operator  R  is  a  con¬ 
traction  on  the  ball  B  while  (136)  says  that  it  maps  the  ball  B  into  itself.  We  now  apply  the 
contraction  mapping  theorem  (Theorem  8  in  Section  2.4)  to  conclude  that  any  sequence  of 
fixed  point  iterates  for  (122)  which  originates  in  B  will  converge  to  a  solution  of  (122).  Since 
is  such  a  sequence,  we  are  done.  □ 


If  ip  is  a  solution  of  (122)  then  the  function  a  defined  by  the  formula 


alx  = 


-»0 

v.  - 


exp (ix£)i/>(£)  d£ 


(137) 


is  clearly  a  solution  of  the  band-limited  integral  equation  (112).  Note  that  because  'ip  e  L1 
the  integral  in  (137)  is  well-defined  and  a  is  an  element  of  the  space  Co  (M)  of  continuous 
functions  which  vanish  at  infinity.  We  record  this  observation  as  follows: 

Theorem  15.  Suppose  that  A  >  0  is  a  real  number.  Suppose  also  that  p  e  L1  (M)  such  that  p  e  L1 
and 


A2 
<  — . 
18 


11  -  To  •  (138) 

Then  there  exists  a  function  a  e  C0  (M)  which  is  a  solution  of  the  integral  equation  (112). 

Remark  2.  Since  a  is  not  necessarily  in  L 1  (M),  the  integral 

exp (-ixp)a(x)  dx  (139) 


(140) 


need  not  exist.  Nor  is  the  existence  of  the  improper  integral 

-ft 


lim 

ft— >00 


ft 


exp(— ixfp)(j{x)  dx 


guaranteed.  However,  when  viewed  as  a  tempered  distribution,  the  Fourier  transform  of  a  exists  and 


18 


is  0;  that  is  to  say, 


rCO  rCO 

j  %l){x)f(x)  dx  =  a(x)f(x)  dx  (141) 

J  —  00  J—  00 

for  all  functions  f  e  S(R). 

In  the  next  section  we  will  prove  that  under  additional  assumptions  on  v,  0  lies  in  L 2  (M).  This  implies 
that  cr  e  L2  (M)  and  ensures  the  convergence  of  the  improper  Riemann  integrals  (140). 


7.  Fourier  estimate 


In  this  section,  we  derive  a  pointwise  estimate  on  the  solution  0  of  Equation  (122)  under 
additional  assumptions  on  the  function  v. 

Lemma  1.  Suppose  that  a  and  C  are  real  numbers  such  that 

0  ^  C  <  a.  (142) 

Suppose  also  that  f  e  L1  (M),  and  that 

1/(01  <  C'exp(— a|^|)  for  all  ^eR.  (143) 

Then 

|exp;[/]K)|«^exp(-a|?|)ii0^expf4^')exp('^:|^  for  all  (  €  K,  (144) 
where  exp^  is  the  operator  defined  in  (119). 


Proof.  Let 


9(0  =  Cexp(-a|0)  (145) 

and  for  each  integer  m  >  0,  denote  by  gm  the  m-folcl  convolution  product  of  the  function  g. 
That  is  to  say  that  g\  is  defined  via  the  formula 

^i(0  =  0(0  (146) 

and  for  each  integer  m  >  0,  gm+i  is  defined  in  terms  of  grn  by  the  formula 

9m+l(0  =  9m  *  g(0-  (147) 

We  observe  that  for  each  integer  m  >  0  and  all  (el, 


9,n({)  =  2  CaC 


2n  J  T(m) 


(148) 


where  Ku  denotes  the  modified  Bessel  function  of  the  second  kind  of  order  v  (see  Section  2.1). 
By  repeatedly  applying  Theorem  1  of  Section  2.1,  we  conclude  that  for  all  integers  m  >  0  and 
all  real  t, 


Km.- 1/2  (£) 


< 


m—  1  / 

Km®  n 

3= 1  V 


2  (J  -  I) 


+  1 


=  K1/2(t)  (  - 


m_1  T  (1±1  +  m  -  1) 

^  r  W) 


(149) 
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We  insert  the  identity 


ki/2  (t)  =  y^exp(-t) 


into  (148)  in  order  to  conclude  that  for  all  integers  m  >  0  and  all  real  numbers  t  >  0, 


^  .  Vvr/tV  ,  r(l±l+m  — 1) 

Km-iM  <  ~y  [2)  exp(-t)  r  - . 

By  combining  (151)  and  (148)  we  conclude  that 


(150) 


(151) 


n  \  m-l  r 

9m(0  <  Cexp(— a|£|)  I  — 


7r  a 


i+ghl 


+  m  —  1 


r(m)r 1 


(152) 


for  all  integers  m  >  0  and  all  £  =/=  0.  Moreover,  the  limit  as  £  — *  0  of  each  side  of  (152)  is  finite 
and  the  two  limits  are  equal,  so  (152)  in  fact  holds  for  all  (et  We  sum  (152)  over  m  =  2,3,... 
in  order  to  conclude  that 


exp^ 


[9]  (0  <  Cexp(-a|£|)  2 


m= 2 


c 


m—  1 


/  l+a,|g| 


+  m  —  1 


W  r(m  +  i)r(m)r 


Cexp(-a|£|)  Yj 


c 


A+alg| 


+  m 


(153) 


for  all  (el.  Now  we  observe  that 

1 


T(m  +  2) 

Inserting  (154)  into  (153)  yields 


1 

€  I  - 


'1  VW  r(m  +  2)r(m  +  i)r 

for  all  m  =  0,1,2, ... . 


exps  [g]  (£)  <  C exp(— a|£|) 


CO  /  C  \  m  r  (  +  m 


m=  1 


2naJ  r(m  +  i)r  ) 


(154) 


(155) 


for  all  £  e  M.  Now  we  apply  the  binomial  theorem  (Theorem  2  of  Section  2.2),  which  is 
justified  since  C  <  a,  to  conclude  that 


exP2  [9]  (0  <  Cexp(— a 


i+“iei 

C  \  “ 


=  C  exp(— a|£|)  exp 


2na 
1  +  a|£| 


(156) 


log 


1  __  c_ 

2tt  a 


for  all  (el.  We  observe  that 


exp(x)  —  1  <  xexp(x)  for  all  x  ^  0, 


and 


1  log 


1  —  X 


^  2x  for  all  0  ^  x  ^  — . 

2tt 


(157) 

(158) 
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By  combining  (157)  and  (158)  with  (156)  we  conclude  that 


exp*  [g]  (£)  <  Cexp(-a|?|)i^!  log 

^  C2  ,th  1  +  a|£| 

<  —  exp(-a|f|) - exp 

Z'K  cl 


1  - 

2na 


exp 


1  +  a  |£| 


log 


1  - 

2tvcl 


c 

27TCL 


“P  ( |l«l 


(159) 


for  all  (ei  Owing  to  (143), 

lexps  [/]  (01  <  exp*  [<)]  (0  for  all  (el.  (160) 

By  combining  this  observation  with  (159),  we  obtain  (144),  which  completes  the  proof.  □ 
Remark  3.  Kummer's  confluent  hypergeometric  function  M(a,  b,  z )  is  defined  by  the  series 


„  r/  ,  .  „  az  (a)2z2  ( a)3z 3 

M^  =  1  +  T  +  W2f  +  \ki  + 


where  (a)n  is  the  Pochhammer  symbol 

T  (a  +  n) 


[  O' IT). 


r(«) 


=  a(a  +  l)(a  +•  2) . . .  (a  +  n  —  1). 


(161) 


(162) 


1  +9a^,2,  —  )  —  1  )  for  all  £  s 
2  7r  a  1  ' 


(163) 


By  comparing  the  definition  of  M{a ,  b ,  z)  with  (153),  we  conclude  that 
|exP;  [/]  «)|«Cexp(— a|£|)(M 

provided 

1/(01  <  c'exp(-a|0)  for  all  (el.  (164) 

The  weaker  bound  (144)  is  sufficient  for  our  immediate  purposes,  but  formula  (163)  might  serve  as  a 
basis  for  improved  estimates  on  solutions  of  Rummer's  equation. 

The  following  lemma  is  a  special  case  of  Formula  (148). 


Lemma  2.  Suppose  that  C  ^  0  and  a  >  0  are  real  numbers,  and  that  f  s  L1  (M)  such  that 

1/(01  <  c  exp  (— a|f|)  for  all  (el.  (165) 

Then 


|/*/(0l  <  C'2exp(-a|^|) 


1  +  a\fi\ 


for  all  £ 


(166) 


We  will  also  make  use  of  the  following  elementary  observation. 


Lemma  3.  Suppose  that  a  >  0  is  a  real  number.  Then 

exp(-a|?l)l?l «  ^5(i) 


for  all  (el. 


(167) 


We  combine  Lemmas  1  and  2  with  (128)  and  (129)  in  order  to  obtain  the  following  key  esti¬ 
mate. 
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Theorem  16.  Suppose  that  r>0,  A>0,  a>0  and  C  ^  0  are  real  numbers  such  that 

0  <  C  <  2aA2. 


Suppose  also  that  f  e  L1  (M)  such  that 

1/(01  <  Cexp(-a|0)  for  all  |f|  <  V2A, 

and  that  v  e  L 1  (M)  such  that 


Kf)l  <  r exp(— a|C|)  forall  £  e 
Suppose  further  that  R  is  the  operator  defined  via  (123).  Then 


I R[f]  (01  <  exp(— a|£|) 
for  all  (eR. 


C 2  /l  +  a|m  f  1  ,  1 


"  +  2^eXP 


C 


47rA  2a 


exp 


C 


47rA2 


10 


+  r 


(168) 

(169) 

(170) 

(171) 


Proof.  We  define  the  operator  via  the  formula 

Rilf]  K)  =  i^[/]*^[/]K) 

and  R-2  by  the  formula 

R2  [f]  (0  =  -4A2exp  *2\Wb[f  ]](£), 

where  Wb  and  114  are  defined  as  in  Section  5.  Then 


(172) 

(173) 


*[/](0  =  *i  [/]  (0  4”  R'2  [f ]  (0+«(0 

for  all  (el.  We  observe  that 

Wb  [/]  (0|  ^  7^  exp(— a|£|)  for  all  (eR. 


By  combining  Lemma  2  with  (175)  we  obtain 

l-Ri  [/]  (01  <  ^exp(-a|^|) 


1  +  a  |£| 


for  all  £  e 


(174) 

(175) 


(176) 


Now  we  observe  that 

\wb [/]  (01  «  L_exp(-a|J|)  for  all  {  6  E.  (177) 

Combining  Lemma  1  with  (177)  yields 

1*2  [/]  (01  «  4^^P(-«kl)  fyffy)  exp  (4^)  <*P  (i^l^l)  <178> 

for  all  £  e  M.  Note  that  (168)  ensures  that  the  hypothesis  (142)  in  Lemma  1  is  satisfied.  We 
combine  (176)  with  (178)  and  (170)  in  order  to  obtain  (171),  and  by  so  doing  we  complete  the 
proof.  □ 


Remark  4.  Note  that  Theorem  16  only  requires  that  /(£)  satisfy  a  bound  on  the  interval  |  -y'2A,  V''2A  | 
and  not  on  the  entire  real  line. 


In  the  next  theorem,  we  use  Theorem  16  to  bound  the  solution  of  (122)  under  an  assumption 
on  the  decay  of  v. 
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Theorem  17.  Suppose  that  A  >  0,  a  >  0  and  T  ^  0  are  real  numbers  such  that 


A  ^  6  max 


Suppose  also  that  v  s  L 1  (M)  such  that 

|i;(£)|  ^  T  exp(— a|£|)  for  all  (el. 
Then  there  exists  a  solution  ofip(0  of  equation  (122)  such  that 


M(OI  <  2r  exp 


for  all  (ei 


(179) 

(180) 

(181) 


Proof  Due  to  (179)  and  (180), 


A2 

MU  ^  — • 

ii  in  -  lg 


(182) 


It  follows  from  Theorem  14  and  (182)  that  a  solution  0(£)  of  (122)  is  obtained  as  the  limit  of 
the  sequence  of  fixed  point  iterates  {VMOI  defined  by  the  formula 

MO  =  v(0  (183) 

and  the  recurrence 


0n+1(O  =  f?[0  n](0-  (184) 

We  now  derive  pointwise  estimates  on  the  iterates  ipn{0  in  order  to  establish  (181). 

Let  {fk\k=o  be  the  sequence  of  real  numbers  be  generated  by  the  recurrence  relation 

&+i  ~2\+T  (185) 

with  the  initial  value 


A)  =  r.  (186) 

It  can  be  established  by  induction  that  (179)  implies  that  this  sequence  is  bounded  above 
by  2r  and  monotonically  increasing,  and  hence  ff.  converges  to  a  real  number  ft  such  that 

0  <  f  <  2T. 


Now  suppose  that  n  ^  0  is  an  integer,  and  that 

hMf)l  <  A»exp(-a|f|)  for  all  f\  <  V2A.  (187) 

When  n  =  0,  this  is  simply  the  assumption  (180).  The  function  bn+i  (£)  is  obtained  from  ipn(0 
via  the  formula 

0n+1(O  =  JR[0](O-  (188) 

We  combine  Theorem  16  with  (188)  and  (187)  to  conclude  that 

|*W.©Nexp(-ak|)  (§  (1^1)  (l  +  M  (jA.)  exp  (AjKl))  +r)  (189) 

for  all  (el.  The  application  of  Theorem  16  is  justified:  the  hypothesis  (168)  is  satisfied  since 

pn  <  2r  ^  2A2a  (190) 
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for  all  integers  n  ^  0.  We  restrict  £  to  the  interval  [— y/2\,  \] 2A]  in  (189)  and  use  the  fact  that 


1  1 
aX  6  ’ 


(191) 


which  is  a  consequence  of  (179),  in  order  to  conclude  that 

L,.  //-m  _ /  f  Pi  ( 1  +  °V2A A  / 1  ,  1  _____  ( 


l^n+i(OI  <  exp(— a|£|)  I 


8  +  2^eXP 


V2A  )  )  +  r 


<exp(-a|®  (xQ  +  V2)(i  +  ^exp( 


Pn 

24ttA 


Pn 

2V2ttA 


for  all  |£|  ^  \/2A.  Now  we  combine  (192)  with  the  inequality 


Pn  <  2T  1 
A  "  A  <  3 


and  the  observation  that 


1  1 


i  +  V2  x  +  ^exp  ^ 


2tt  1  \72n 


(192) 


(193) 


(194) 


in  order  to  conclude  that 


IVW,(«)I  «  eXp(-0|{|)  Q  +  V2)  Q  +  4  exp  (W)  exp  (; 


«  (  ^  +  rJexp(-a|{|) 


(195) 


=  Pn+i  exp(— a|£|) 

for  all  |£|  ^  V2A.  We  conclude  by  induction  that  (187)  holds  for  all  integers  n  ^  0. 

The  sequence  {pn{P)}  converges  to  h(£)  in  L1  (M)  norm  (and  hence  a  subsequence  of  ■0n(£) 
converges  to  £;(£)  pointwise  almost  everywhere)  and  (187)  holds  for  all  integers  n  ^  0.  More¬ 
over,  for  all  integers  n  ^  0,  pn+  \  ^  2 T.  From  these  observations  we  conclude  that 

hHOI  <  2r  exp(-a|£|)  (196) 

for  almost  all  |£|  ^  \/2A.  We  also  observe  that  ip(£)  is  a  fixed  point  of  the  operator  R,  so  that 

m  =  RW(o  (i  97) 

for  all  (eK.  Clearly, 

R[f](0  =  R[g](0  (198) 

for  all  £  e  R.  if  /(£)  =  <?(£)  for  almost  all  |£|  ^  y/2X,  so  (196)  in  fact  holds  for  all  |£|  ^  ^/2X. 

We  now  apply  Theorem  16  to  the  function  ^(£)  (which  is  justified  since  2 T  <  2A2a)  to  conclude 
that 


1^(01  <  exp(— a|£|)  (  (  - - 


4T2  /l  +  a|£|\  (1  ,  1 


8  +  2^eXP  4^ 


l£l  +  r  (199) 


for  all  (el.  Note  the  distinction  between  (187)  and  (199)  is  that  the  former  only  holds  for  all 
£  in  the  interval  |  -y2A,  V''2A  |,  while  the  later  holds  for  all  £  on  the  real  line.  It  follows  from 
(179)  that 


11  T  1 

Aa  <  6  and  A  <  6 


(200) 
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We  insert  these  bounds  into  (199)  in  order  to  conclude  that 

«?)|  «  rexpMfl)  (T  (i^)  (5  +  ^exp  (W)  exp  (j^l))  +  l)  (201) 
for  all  (el.  Now  we  observe  that 


2vr 


exp 


72vr 


1 

<  6’ 


which,  when  combined  with  (201),  yields 


1^(01  <  r exp(-a^l) 


1  2j£| 

9  3A 


1  1 
8  +  6eXP 


12ttA  1 


|f|  +  1  for  all  £  e 


By  rearranging  the  right-hand  side  of  (203)  as 

/  1  1  /  1 
r  exp(— a|£|)  —  +  —  exp 


72  56 


12ttA  1 


m  + 


lei  ^  lei 

12A  +  9A  6XP 


12ttA  1 


lei  +i 


=  Texp  (  -  (  a  —  -  )  |£|  )  exp  (  -T|f| 


A 


A1 


- 1 - exp 

72  56  1 


12vrA 


lei 


lei  ^  lei 

+  l2A  +  9AeXP 


and  applying  Lemma  3,  we  arrive  at  the  inequality 

MO  I  <  rexP  (~  (a  -  t)  \t\)  (ik  +  h  + 


12ttA  1 


+ 


10  +1 


72  56  12exp(l)  33exp(l) 


(202) 


(203) 


+  1  for  all  £  e 


from  which  (181)  follows  immediately 


(204) 

□ 


Suppose  that  f  e  L1  (M)  is  a  solution  of  (122).  Then  the  function  a  defined  by  the  formula 

rCO 

a(x)  =  exp (ix£)i/;(€)  (205) 

J  —  00 

is  a  solution  of  the  integral  equation  (112).  However,  the  Fourier  transform  of  (205)  might 
only  be  defined  in  the  sense  of  tempered  distributions  and  not  as  a  Lebesgue  or  improper 
Riemann  integral.  If,  however,  we  assume  the  function  p  appearing  in  (122)  is  an  element 
of  L 1  (M)  and  impose  the  hypotheses  of  Theorem  18  on  the  Fourier  transform  of  p,  then  u>  e 
L 2  (M),  from  which  we  conclude  that  a  is  also  an  element  of  L2  (M).  In  this  event,  there  is 
no  difficulty  in  defining  the  Fourier  transform  of  a.  We  record  these  observations  in  the 
following  theorem. 


Theorem  18.  Suppose  that  there  exist  real  numbers  A  >  0,  T  >  0  and  a  >  0  such  that 

A  >  6  max  |r,  -  j  .  (206) 

Suppose  also  that  p  e  L1  (M)  such  that  and 

|p(OI  ^rexp(-alel)  foralUeR.  (207) 

Then  there  exists  a  solution  a(x)  of  the  integral  equation  (112)  such  that 

\£\)  for  all  (eK  (208) 


|a(0|  <  2 T  exp 


1 
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Remark  5.  The  bound  (208)  implies  that  a  decays  faster  than  any  polynomial,  from  which  we  conclude 
that  a  is  infinitely  differentiable. 


8.  Approximate  solution  of  the  original  equation 


We  would  like  to  insert  the  solution  a  of  (112)  into  the  original  equation  (74).  However,  we 
have  no  guarantee  that  a  is  in  L 1  (M),  nor  do  we  expect  that  a(+2A)  =  0.  As  a  consequence, 
the  integrals  defining  T  [cr]  might  not  exist. 

To  remedy  this  problem,  we  define  a  "band-limited"  version  ab  of  a  by  the  formula 

*b(Z)  =  Z(ZM),  (209) 


where  &(£)  is  the  function  used  to  define  the  operator  Tb.  We  observe  that  there  is  no  difficulty 
in  applying  T  to  ab  since 


cp,(+ 2A)  =  0. 

(210) 

Moreover,  Tb[a ]  =  T  [cr6],  so  that 

a(x)  =  S[T  [crb]]  ( x )  +  p(x) 

for  all 

(211) 

Rearranging  (211),  we  obtain 

(Tb(x)  =  S[T  [cr6]  ]  (x)  +  pb(x) 

for  all 

(212) 

where  pb[x)  is  defined  the  formula 

Pb{x)  =  p(x)  +  ab(x) 

-a(x). 

(213) 

Using  (208)  and  (213),  we  conclude  that  under  the  hypotheses  of  Theorem  18, 

\\p-PbU  T 
< 


< 


< 


Wb-cr  ||i 

2T  exp 
J|€I>A 
4T 


A 


lei  df 


a  A 


r  exp  (  -  (  a  -  -  ]  A 


(214) 


24T 
5  a 


exp 


pa\ 

o 


Together  Theorem  18  and  (214)  imply  Theorem  12. 


9.  Backwards  error  estimate 

In  this  section,  we  prove  Theorem  13. 

Although  both  p  and  pb  are  defined  on  the  real  line,  we  are  only  concerned  with  solutions  of 
(1)  on  the  interval  [0, 1],  and  so  we  only  require  estimates  there.  Accordingly,  throughout  this 
section  we  use  ||  ■  ||oo  to  denote  the  ([0, 1])  norm. 


26 


We  omit  the  proof  of  the  following  lemma,  which  is  somewhat  long  and  technical  but  entirely 
elementary  (it  can  be  established  with  the  techniques  found  in  ordinary  differential  equation 
textbooks;  see,  for  instance.  Chapter  1  of  [4]). 


Lemma  4.  Suppose  that  q  :  M.  ->  M  is  infinitely  differentiable  and  strictly  positive,  that  p(t )  is  defined 
by  (56),  and  that  there  exist  real  numbers  yi  >  0  and  q2  >  0  such  that 

r/i  ^  q(t )  ^  y2  for  all  0  ^  t  ^  1,  (215) 

and 


Let 


\p(t)\ ,  k'COI  <  m  for  all  0  ^  t  <  1. 


k  =  2D[  —  j  +8^  +  10— +  1 

\mj  vi 


and  suppose  also  that  e  >  0  is  a  real  number  such  that 


^  r6 
6  <  2  ‘ 


Suppose  furthermore  that  pb  :  [0, 1]  ->  M  is  an  infinitely  differential  function  such  that 

\\p~Pb\\co  c  eexp (-k). 

Then  there  exists  an  infinitely  differentiable  function  qb  :  [0, 1]  — >•  R.  such  that 


Qb(t)  l  4  \qb(t)J  qb(t) 


(216) 

(217) 

(218) 

(219) 

(220) 


and 


q  -  g&IU  <  e- 


(221) 


We  derive  a  bound  on  the  change  in  solutions  of  the  ordinary  differential  equation  (1)  when 
the  coefficient  q(t)  is  perturbed. 

Lemma  5.  Suppose  that  A  >  0,  e  >  0,  ry  >  0  and  y2  >  0  are  real  numbers.  Suppose  also  that 
q  :  [0, 1]  ->  M  is  a  continuously  differentiable  function  such  that 

rji  ^  q(t)  ^  y2  for  all  0  <  £  <  1,  (222) 

that  p(t)  is  defined  by  the  formula  (56),  and  that 

\p(t)\  ^  r/2  for  all  0  <  £  <  1.  (223) 

Suppose  furthermore  that  qb  :  [0, 1]  — »  M  is  a  continuously  differentiable  function  such  that 

\q{t)  —  Qb(t) I  <  2Xexp  (  )  e  foral1  0  ^  ^  L  (224) 

If  z(t)  is  a  solution  of  the  ordinary  differential  equation 

z"(t)  +  A 2q(t)z(t)  =  0  for  all  0  <  f  ^  1  (225) 

and  z0(t )  is  the  unique  solution  of  the  ordinary  differential  equation 

z'lft)  +  X2qb(t)z0(t )  =  0  for  all  0  ^  t  ^  1  (226) 
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such  that  z0(0)  =  z(0)  and  z'0( 0)  =  z'(0),  then 

\\z  -  Zq || oo  ^  e|k||oo- 


(227) 


Proof.  We  start  by  observing  that  the  function  'ip(t)  =  z0(t )  —  z(t)  is  the  unique  solution  of  the 
initial  value  problem 

ip"(t)  +  A  2q(t)ip(t)  =  A2/(t)  for  all  0  ^  t  ^  1 

m  =  v>'(o)  =  o, 


(228) 


where 


f(t)  =  (q(t)  -  qo(t))  (f(t)  +  z{t)) . 


(229) 


We  now  apply  the  well-known  Liouville-Green  transformation  to  (228)  in  two  steps.  First, 
we  introduce  the  function 

4>{f)  =  (<?O))1/4,0W,  (230) 

which  is  the  solution  of  the  initial  value  problem 

<!>"(*)  -  +  A 2q{t)<t>{t)  =  (g(f))1/4  A 2 f{t)  -  ^ q(t)p(t)(f>(t )  for  all  0  <  t  <  1 

HO)  =  0'(O)  =  0. 

Next  we  introduce  the  change  of  variables 

Vq(u)  du, 

Jo 


x{t)  = 


(232) 


which  transforms  (231)  into 

(x)  +  \2<p(x 

m  =  0'(o)  =  o, 

where 


4>"(x)  +  A 2f(x)  =  (q(x))  3/4  A 2f(x)  —  -p(x)f(x)  for  all  0  ^  x  ^  x\ 

JL  ) 


xi  =  V q(u)  du. 


(234) 

We  now  use  the  Green's  function  for  the  initial  value  problem  (233)  obtained  from  the  Liouville- 
Green  transform  to  conclude  that  for  all  0  ^  x  ^  x\, 


4>{x)  = 


Jo 


sin(A(x  —  y)) 
A 


(q(y))  3/4  A 2f(y )  -  - p(y)<t>(y)  )  dy. 


(235) 


By  inserting  (229)  into  (235)  and  we  obtain  the  inequality 


A  772 


Vi 


A 


I  <  Ik  -  qb\\ao—  +  t  \<f>(y)\  dy  +  ^n\\q  -  qb\\ao\\z\\ao  for  all  0  ^  x  <  xi.  (236) 


3/4  I 

Vi 


We  apply  Gronwall's  inequality  (Theorem  9  in  Section  2.5)  to  (236)  in  order  to  conclude  that 

A 


3/4 

Vi 


\q  -  gblloolklloo  exp  I  (  Ik  -  g&||oo—  +  ~r  )  x  )  for  all  0  <  x  ^  xx . 

T)i  4 


(237) 
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(238) 


Combing  (230),  (237)  and  the  observation  that 


X\  =  \f  q(u)  du  ^ 


yields  the  inequality 


A 


Vi 


A 


rh 


\<l>(t)\  <  —h  -  Qb  II  oo  exp  —  \\q  -  qb\\ooV22  exP 


1/2 


Vi 


3/4' 

% 


for  all  0  ^  t  ^  1.  Now 


A 


Vi 


A 


—  Ik  -  ^llooexp  —\\q  -  qb\\oDV2/2  exP 


Vi 


3/4' 

V2 


<  e 


if  and  only  if 


—vl/2\ \q  ^  Qb\\oo  <  W0  (  e^exp 

V  i 


1/2, 


3/4' 

V2 


(239) 


(240) 


(241) 


where  W0  is  the  branch  of  the  Lambert  W  function  which  is  real-valued  and  greater  than  —1 
on  the  interval  [—1/e,  oo)  (see  Section  2.3).  According  to  Theorem  5, 


A 


-d22h  -  <?b||oo  <  k  I  e??2/Z  eXP  (  ^ 

implies  (241).  We  algebraically  simplify  (242)  in  order  to  conclude  that 


!/2, 


3/4' 

V2 


ii  ii  ^  1 

Ik  -  9b||oo  <  2Xexp  I  “ 

implies  (239). 

By  combining  Lemmas  4  and  5  we  obtain  Theorem  13. 


3/4' 

V2 


(242) 

(243) 

□ 


10.  Numerical  experiments 

In  this  section,  we  describe  numerical  experiments  which,  inter  alia,  illustrate  one  of  the  im¬ 
portant  consequences  of  the  existence  of  nonoscillatory  phase  functions.  Namely,  that  a  large 
class  of  special  functions  can  be  evaluated  to  high  accuracy  using  a  number  of  operations 
which  does  not  grow  with  order. 

Although  the  proof  of  Theorem  12  suggests  a  numerical  procedure  for  the  construction  of 
nonoscillatory  phase  functions,  we  utilize  a  different  procedure  here.  It  has  the  advantage 
that  the  coefficient  q  in  the  ordinary  differential  equation  (1)  need  not  be  extended  outside  of 
the  interval  on  which  the  nonoscillatory  phase  function  is  constructed.  A  paper  describing 
this  work  is  in  preparation. 

The  code  we  used  for  these  calculations  was  written  in  Fortran  and  compiled  with  the  Intel 
Fortran  Compiler  version  12.1.3.  All  calculations  were  carried  out  in  double  precision  arith¬ 
metic  on  a  desktop  computer  equipped  with  an  Intel  Xeon  X5690  CPU  running  at  3.47  GFIz. 
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10.1.  A  nonoscillatory  solution  of  the  logarithm  form  of  Kummer's  equation. 

In  this  experiment,  we  illustrate  Theorem  12  in  Section  4.  We  first  construct  a  nonoscillatory 
solution  r  of  the  logarithm  form  of  Kummer's  equation 

r"{t)  —  ^  fit))2  +  4A2  (exp  (r(t))  —  q(t ))  =  0  (244) 

on  the  interval  [—1, 1],  where  A  =  1,000  and  q  is  the  function  [—1, 1 1  — >  M  defined  by  the 
formula 

q(t)  =  ^3  +  -  +11Q^2  +  f3  cos(5f)^  .  (245) 

Then  we  compute  the  500  leading  Chebyshev  coefficients  of  q  and  r. 

We  display  the  results  of  this  experiment  in  Figures  1  and  2.  Figure  1  contains  plots  of  the 
functions  q  and  r,  while  Figure  2  contains  a  plot  of  the  base-10  logarithms  of  the  absolute 
values  of  the  leading  Chebyshev  coefficients  of  q  and  r. 

We  observe  that,  consistent  with  Theorem  12,  the  Chebyshev  coefficients  of  both  r  and  q 
decay  exponentially,  although  those  of  r  decay  at  a  slightly  slower  rate. 


10.2.  Evaluation  of  Legendre  polynomials. 


In  this  experiment,  we  compare  the  cost  of  evaluating  Legendre  polynomials  of  large  order 
using  the  standard  recurrence  relation  with  the  cost  of  doing  so  with  a  nonoscillatory  phase 
function. 


For  any  integer  n  ^  0,  the  Legendre  polynomial  Pn{x)  of  order  n  is  a  solution  of  the  second 
order  differential  equation 

(1  -  t2)y"(t)  -  2 ty'(t)  +  n(n  +  1  )y(t)  =  0.  (246) 


Equation  (246)  can  be  put  into  the  standard  form 

,h(4.\  ,  f±  +  n-nt2 -n2(t2 -1) 

*  {t)  +  ( - qrw - 

by  introducing  the  transformation 


fit)  =  o 


(247) 


'fit)  =  Vl  -  t2  y(t). 


(248) 


Legendre  polynomials  satisfy  the  well-known  three  term  recurrence  relation 

(n  +  1  )Pn+i(t)  =  (2 n  +  1  )tPn(t)  -  nPn-i(t).  (249) 

See,  for  instance,  [11]  for  a  discussion  of  the  these  and  other  properties  of  Legendre  polyno¬ 
mials. 


For  each  of  9  values  of  n,  we  proceed  as  follows.  We  sample  1000  random  points 

ti,  t2,  ■  ■  ■  ,  fiooo  (250) 

from  the  uniform  distribution  on  the  interval  (—1, 1).  Then  we  evaluate  the  Legendre  polyno¬ 
mial  of  order  n  using  the  recurrence  relation  (249)  at  each  of  the  points  C,  f2,  •  •  • ,  fiooo-  Next, 
we  construct  a  nonoscillatory  phase  function  for  the  ordinary  differential  equation  (247)  and 
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use  it  evaluate  the  Legendre  polynomial  of  order  n  at  each  of  the  points  ti,  t2, . . . ,  tiooo-  Fi¬ 
nally,  for  each  integer  j  =  1, ... ,  1000,  we  compute  the  error  in  the  approximation  of  Pn(tj ) 
obtained  from  the  nonoscillatory  phase  function  by  comparing  it  to  the  value  obtained  us¬ 
ing  the  recurrence  relation  (we  regard  the  recurrence  relation  as  giving  the  more  accurate 
approximation). 

The  results  of  this  experiment  are  shown  in  Table  1.  There,  each  row  correponds  to  value  of 
n.  That  value  is  listed  n,  as  is  the  time  required  to  compute  each  phase  function  for  that  value 
of  n,  the  average  time  required  to  evaluate  the  Legendre  polynomial  of  order  n  using  the 
recurrence  relation,  the  average  cost  of  evaluating  the  Legendre  polynomial  of  order  n  with 
the  nonoscillatory  phase  function,  and  the  largest  of  the  absolute  errors  in  the  approximations 
of  the  quantities 

Pn(t l),  Pn(t2)r  •  •  •  >  -Pn(t10oo) 

obtained  via  the  phase  function  method. 

This  experiment  reveals  that,  as  expected,  the  cost  of  evaluating  Pn(t)  using  the  recurrence 
relation  (249)  grows  as  0(n )  while  the  cost  of  doing  so  with  nonoscillatory  phase  function  is 
independet  of  order. 

However,  it  also  exposes  a  limitation  of  phase  functions.  The  values  of  Pn(t )  are  obtained  in 
part  by  evaluating  sine  and  cosine  of  a  phase  function  whose  magnitude  is  on  the  order  of  n. 
This  imposes  limitations  on  the  accuracy  of  the  method  due  to  the  well-known  difficulties  in 
evaluating  periodic  functions  of  large  arguments. 

Figure  3  contains  a  plot  of  the  nonoscillatory  phase  function  for  the  equation  (247)  when 

n  =1,000,000. 


10.3.  Evaluation  of  Bessel  functions. 

In  this  experiment,  we  compare  the  cost  of  evaluating  Bessel  functions  of  integer  order  via 
the  standard  recurrence  relation  with  that  of  doing  so  using  a  nonoscillatory  phase  function. 

We  will  denote  by  Jv  the  Bessel  function  of  the  first  kind  of  order  u.  It  is  a  solution  of  the 
second  order  differential  equation 

t2y"(t)  +  ty'(t)  +  (t2-v2)y(t)  =  0,  (251) 

which  can  be  brought  into  the  standard  form 

V(t)  +(l-  -  ~ 1/4 )  m  =  o  (252) 

via  the  transformation 

ip(t)  =  Vty(t).  (253) 

An  inspection  of  (252)  reveals  that  Ju  is  nonoscillatory  on  the  interval 

(o,  ^  V4/A  -  l)  (254) 
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and  oscillatory  on  the  interval 


^  V4z/2  -  1,  oo^  .  (255) 

The  Bessel  functions  satisfy  the  three-term  recurrence  relation 

Jv+i(t)  =  y4(t)-4-i(l).  (256) 

The  recurrence  (256)  is  numerically  unstable  in  the  forward  direction;  however,  when  evalu¬ 
ated  in  the  direction  of  decreasing  index,  it  yields  a  stable  mechanism  for  evaluating  Bessel 
functions  of  integer  order  (see,  for  instance.  Chapter  3  of  [11]). 

For  each  of  9  values  of  n,  we  proceed  as  follows.  First,  we  sample  1000  random  points 

ti,  t2,  ■  ■  ■ ,  fiooo  (257) 

from  the  uniform  distribution  on  the  interval  [2n,3n].  We  then  use  the  recurrence  relation 
(256)  to  evaluate  the  Bessel  function  Jn  of  order  n  at  the  points  t\,t2,  ■  . .  , tiooo-  Next,  we 
construct  a  nonoscillatory  phase  function  for  the  equation  (253)  on  the  interval  [2 n,  3n]  and 
use  it  to  evaluate  Jn  at  the  points  A,  t2, . . . ,  t1000.  Finally,  for  each  integer  j  =  1, . . . ,  1000, 
we  compute  the  error  in  the  approximation  of  Jn(tj )  obtained  from  the  nonoscillatory  phase 
function  by  comparing  it  to  the  value  obtained  using  the  recurrence  relation  (once  again  we 
regard  the  recurrence  relation  as  giving  the  more  accurate  approximation). 

The  results  of  this  experiment  are  displayed  in  Table  2.  There,  each  row  corresponds  to  one 
value  of  n.  In  addition  to  that  value  of  n,  it  lists  the  time  required  to  compute  the  phase  func¬ 
tion  at  order  n,  the  average  cost  of  evaluating  .Jn  using  the  recurrence  relation,  the  average 
cost  of  evaluating  it  with  the  nonoscillatory  phase  function,  and  the  largest  of  the  absolute 
errors  in  the  approximations  of  the  quantities 

Jn\p2)i  •  •  ■  i  Jn{t lOOo) 

obtained  via  the  phase  function  method. 

We  observe  that  while  the  cost  of  evaluating  Jn  using  the  recurrence  relation  (256)  grows  as 
0{n),  the  time  taken  by  the  nonoscillatory  phase  function  approach  scales  as  0(1).  We  also 
note  that,  as  in  the  case  of  Legendre  polynomials,  there  is  some  loss  of  accuracy  with  the 
phase  function  method  due  to  the  difficulties  of  evaluating  trigonometric  functions  of  large 
arguments. 


11.  Conclusions 

We  have  shown  that  the  solutions  of  a  large  class  of  second  order  differential  equations  can 
be  accurately  represented  using  nonoscillatory  phase  functions. 

We  have  also  presented  the  results  of  numerical  experiments  which  demonstrate  one  of  the 
applications  of  nonoscillatory  phase  functions:  the  evaluation  of  special  functions  at  a  cost 
which  is  independent  of  order.  An  efficient  algorithm  for  the  evaluation  of  highly  oscillatory 
special  functions  will  be  reported  at  a  later  date. 

A  number  of  open  issues  and  questions  related  to  this  work  remain.  Most  obviously,  a  fur¬ 
ther  investigation  of  the  integral  equation  (74)  and  the  conditions  under  which  it  admits  an 
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exact  solution  is  warranted.  Moreover,  there  are  applications  of  nonoscillatory  phase  func¬ 
tions  beyond  the  evaluation  of  special  functions  which  should  be  explored.  And,  of  course, 
the  generalization  of  these  results  to  higher  dimensions  is  of  great  interest.  The  authors  are 
vigorously  pursuing  these  avenues  of  research. 
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Figure  1:  The  function  q  defined  by  formula  (245)  in  Section  10.1  (solid  line)  and  the  corre¬ 
sponding  solution  r  of  the  logarithm  form  of  Rummer's  equation  (53)  when  A  =  1,000  (dotted 
line). 
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Figure  2:  The  base-10  logarithms  of  the  leading  Chebyshev  coefficients  of  the  function  q  de¬ 
fined  by  formula  (245)  in  Section  10.1  (solid  line)  and  of  the  associated  nonoscillatory  solution 
r  of  equation  the  logarithm  form  of  Rummer's  equation  (53)  when  A  =  1,000  (dotted  line). 
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n 

Phase  function 
construction  time 

Avg.  phase  function 
evaluation  time 

Avg.  recurrence 
evaluation  time 

Largest 
absolute  error 

101 

1.55x10- 

-1  secs 

1.29  x  10~6  secs 

5.82  xlO-8  secs 

5.16x  10~14 

102 

1.76x10- 

-1  secs 

1.29  x  10~6  secs 

9.73  x  10 w  secs 

1.59  xl0~13 

o 

CO 

1.57x10- 

-1  secs 

1.29  x  10”6  secs 

1.03  xl0~5  secs 

6.13x  10"13 

104 

1.55x10- 

1  secs 

1.29  x  10~6  secs 

1.04  x  10-4  secs 

1.20x  10~12 

105 

1.56x  10- 

1  secs 

1.31  x  10~6  secs 

1.04  x  10-3  secs 

9.79  xhT12 

106 

1.58x10- 

-1  secs 

1.40  x  10”6  secs 

9.81  xl0~3  secs 

2.40  xlO"11 

107 

1.65x  10- 

1  secs 

1.40  x  10~6  secs 

9.69  xl0~2  secs 

8.59xl0-11 

108 

1.87x10- 

-1  secs 

1.42  x  10~6  secs 

9.68  xlO-1  secs 

1.71  xlO”10 

109 

2.05x10- 

-1  secs 

1.34  x  10"6  secs 

9.68  xlO-0  secs 

6.11  xlO”10 

Table  1:  The  evaluation  of  Legendre  polynomials.  A  comparison  of  the  time  required  to 
evaluate  the  Legendre  polynomial  of  order  n  using  the  standard  recurrence  relation  and  the 
time  necessary  to  evaluate  it  using  a  nonoscillatory  phase  function.  The  recurrence  relation 
approach  scales  as  0(n)  while  the  phase  function  approach  scales  as  0(1). 


Figure  3:  A  phase  function  for  Legendre's  differential  equation.  A  plot  of  the  nonoscillatory 
phase  function  associated  with  Legendre's  equation  (246)  at  order  n  =  1,000,000.  It  is  suffi¬ 
cient  to  construct  the  phase  function  on  the  interval  [0, 1)  due  to  the  symmetry  properties  of 
Legendre's  differential  equation. 
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n 

Phase  function 
construction  time 

Avg.  phase  function 
evaluation  time 

Avg.  recurrence 
evaluation  time 

Targest 
absolute  error 

101 

5.23x10- 

~4  secs 

1.30  xlO”6  secs 

1.99  x  10~6  secs 

2.81  xlO”14 

102 

5.39x10 

~4  secs 

1.31  x  10”6  secs 

7.29  x  10~6  secs 

7.85  xlO”14 

103 

5.36x10 

~4  secs 

1.37xl0”6  secs 

4.87x  10~5  secs 

2.40  xl0“13 

104 

5.52x10- 

~4  secs 

1.33 xlO”6  secs 

4.35  xl0~4  secs 

1.01  xlO”12 

105 

5.46  x  10 

~4  secs 

1.49  x  10”6  secs 

4.11  x  10-3  secs 

3.18xl0”12 

106 

5.81x10 

~4  secs 

1.44  x  10"6  secs 

4.24 xlO-2  secs 

8.57xl0~12 

107 

6.41x10 

~4  secs 

1.45  xl0~6  secs 

4.36  xlO-1  secs 

5.98xl0”41 

108 

7.00x10 

~4  secs 

1.35 xlO”6  secs 

4.39  x  10+°  secs 

1.14xl0”10 

109 

1.26  x  10+°  secs 

1.41  x  10~6  secs 

4.42  x  10+1  secs 

2.43  xlO“10 

Table  2:  The  evaluation  of  Bessel  functions.  A  comparison  of  the  time  required  to  evaluate 
the  Bessel  function  Jn  using  the  standard  recurrence  relation  with  that  required  to  evaluate 
it  using  a  nonoscillatory  phase  function.  All  of  the  points  at  which  Jn  was  evaluated  were  in 
the  interval  [2 n,  3n] .  The  recurrence  relation  approach  scales  as  0(n)  in  the  order  n  while  the 
time  required  by  the  phase  function  method  is  0(1). 
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